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a m n non 


LINER is a system of Fortran 77 codes which perforins a 2D analysis of acoustic wa ^ 
•opagation and noise suppression in a rectangular channel with a contrnuons ner 
le top wall. This new implementation is designed to streamline the usage of t e severa 
Ls mahmg up LINER, resulting in a useful design too,. Ma.r input parame „ 
w in two main data hies, inline and num.pnn. Output data appear m th 
,f ASCII files as well as a choice of GNUPLOT graphs. Section 2 briefly deserrbes 
JZ model, sectron 3 discusses the numerical methods, Section 4 gives a detarled 
recount of program usage, including input formats and graphical options. A sample run 
„ abo provided. Finally, Section 5 briefly describes the individual program es. 

*2 Physical Model 

The physical model is based on the 2D compressible, inviscid fluid flow m a rec an- 
gular channel, as depicted in Fig. 1. We employ Cartesian coordinates 
„ the streamwise direction and y the normal direction. The governing e q uat.ons are 

2D Euler equations, 

Continuity : + V ’ = 0 

. . . j o Pirinffen Report and Associated Software 

Copyright 1999 by R. S. Reichert and S. Binngen, Kep 
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Momentum : 


p V- + p\ • VV = — Vp 

dt 


Energy : + V • EY — —V • (pV) 

where V = (u, v) T is the velocity, p is density, p is the pressure, and E is the total energy. 
These are numerically integrated, along with the equation of state, 

E = JL r j + V 2 

where 7 is the ratio of specific heats, which we shall take as 7 = 1-4. Nonreflective inflow 
and outflow boundary conditions are imposed. Boundary conditions at the normal walls 
(y _ ±y max ) are rigid slip, allowing nonzero tangential velocity while maintaining zero 
normal velocity. The initial conditions consist of a quiescent field, with an optional 
y-directed bias flow superimposed. Although this model is inviscid, a small artificial 

viscosity is used to stabilize the difference scheme. 

The liner is incorporated as a source term in the momentum equations, rather than 
through boundary conditions. This is done semi-empirically, with impedance parameters 
chosen to match experimental data on the liner materials. As shown in Fig. 2, a steady 
bias flow through the liner is also modelled. See Ref. 1-2 for a detailed description 
of finer characteristics. Physically, the finer is composed of a face sheet and septum, 
sandwiching a honeycomb core. A constant bias velocity may be directed through the 

finer, if desired. 

3. Numerical Method 

The Euler equations are solved in space and time using the explicit (2,4) scheme of 
Gottlieb & Turkel [3], with artificial viscosity incorporated as a sixth-order source term. 
Spatial boundary conditions are implemented using the method of Thompson [4]. 


2 


4. Computer Programs 


4.1 General Description 

LINER consists of a suite of F77 codes, which may be classified functionally as data 
input (relatively permanent or relatively ephemeral), computational, and data output 
(numerical or graphical). An entire program execution is supervised by a single Unix 

script, called RUN. 

4.2 Input 

Input parameters are fouad in input.inc and nam.prm. The liner geometry is gen- 
erated in Inrgeom.f, but the user need not edit this file. See section 4.5 (Sample Run) 
for the specific contents of the two input files. 

4.3 Output 

Sound pressure level is written to ASCII file sound. idl, while power per span is 
written to file powspan. 

4.4 Graphical Options 

AU plotting is carried out in the GNUPLOT system, with a single Gnuplot macro 
processing a number of ASCII datafiles. Plots are automatically made (using GHOSTVIEW/, 
of Transmitted Power per Span and Sound Pressure Level (SPL) y-profiles at three se- 
lected n- values. To page through plots, type in “q" after each plot. The four plots are 
written to files pl.ps, p2.ps, p3.ps, and p4.ps. The user may produce additional SPL 
profiles by running program splx./with xl, x2, and x3 set as desired in input.inc. 
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4.5 Sample Run 


As a specific case, let us consider a 61 cm x 19.25 cm domain, discretized by a 
489 x 155 grid in x and y, Le. nx = 487, ny = 153. With these choices, dx = dy = 
0.125 cm. The lower finer sheet is located at y index il = 122, with upper sheet (septum) 
at iu = 142, corresponding to a total finer depth of 4 cm. The bias velocity, scaled by 
the sound speed, is negative blowing out of the finer. Here we take i/ 6ias = -0.01. The 
equations are advanced a total of nmx = 6000 steps with timestep dt = 0.02. The input 
disturbance is a plane wave of nondimensional frequency u = 0.23099946 (/ = 1250 Hz) 
and normalized amplitude A = 0.002 (14MB). Specifically, / is scaled by u r /l r , where 
u r is sea level sound speed and l r = 1 cm; A is scaled by p r u 2 r , where p r is sea level 
density. Parameters nmx, dx, and dt are to be specified by the user in file num.prm. 

The finer impedance parameters (defined in Ref. 1) are specified by the user in 
file input.inc. This file also includes parameters dy, nx, ny, il, tu and v bzas . For 
this example, the finer impedance parameters were chosen to be al = 0.0283, au = 
0.0905, bl = 21.8, bu = 233.0, xclsp = -1.14, and xcusp = -1.11338. 

The Sound Pressure Level is to be plotted vs y at xl = 8.5 cm, x2 = 17.0 cm, 
and x3 = 52.5 cm, as set in input.inc. SPL profiles may easily be generated for any 
other desired values of x without rerunning the whole code, by resetting xl,x2, and x3 
in input.inc, then recompiling and executing splx.f. Figure 3 shows the resulting plot of 
Power Drop, while Fig. 4 presents the three SPL plots. 

5. Glossary of Computer Programs 

ff.f: Creates forcing function field to prevent the mean bias field from evolving in time. 
Writes file ff. in. 

goplot: GNUPLOT macro to make plots of power drop and SPL. 
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input-inc: This is the basic input file, to be edited by the user. It is included in all 
FORTRAN codes. 

linchn-orig.f: This is the workhorse of the system, which numerically solves the differ- 
ence equations representing the flow. 

linchn_ppav.f: Same as It nchn.orig.f, except that time-averaging is applied to the fluc- 
tuating (acoustic) field. Writes final output files. 

lnrgeom.f: Constructs liner forcing function distribution throughout domain, 
mnbias.f: Creates solution field for any desired mean flow, such as a bias flow. Writes 
file rstrt.in. 

num.prm: Contains basic run parameters, to be edited by the user. Warning:, do not 
add any additional lines of comments to this file. 

run: Unix script which supervises an entire run. 

spl.f: Postprocessor to generate plot files for goplot 


References 

1. Reichert, R. S.,Ph. D. Thesis, University of Colorado Dept, of Aerospace Engineering 

Sciences, 1998. 

2. Reichert, R. S. and Biringen, S., AIAA 97-1650 (1997). 

3. Gottleib, D. and Turkel, E. Math. Comp. 115, 43 (1976). 

4. Thompson, K. W., J. Comp. Phys. 68, 1 (1987). 
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Nonreflective 

Outflow 













SPL (dB) 





Script to run L 
Generate y-grid 
install test fo 
n newstrt.bin 
grid . out 
7 -u grid.f -o 


o 'finished running 
linchn_orig . f -02 
7 line . f -o gin 

o ' finished running 
newstrt.bin rstrt . i 
pp.av jetdat.bin je 
nura . prra numsave . prm 
o 'starting addn' 
addn.f -o ga 


ho 'finished gl 
newstrt.bin je 
7 linchn_ppav . f 
v 

ho 'finished ru 
num.prm numxx. 
■ numsave . prm n\. 
7 splx.f -o gg 




input . inc 

Main input file, to be edited by the user, 
grid constants (to be changed by user) 

nx = number of interior x gridlines (nx + 2 total, including boundaries) 
ny = number of interior y gridlines (ny+2 total, including boundar 
il = y index of lower porous sheet (face sheet) 
iu = y index of upper porous sheet (septum) 
dely = dx = dy = mesh spacing 

integer nx,ny,il,iu 

parameter (nx=4 87, ny=15 3, dely = . 125 , il=122 , iu=142) 

physical constants (not to be changed) 
real gam, garni 

parameter (gam=l . 4 , garni = 0.4) 

rba and pba are quiescent medium nondim density and pressure 

real rba, pba ___. 

parameter (rba=l . 0 ,pba=0 . 714285714286) 

Print interval 
integer npr 
parameter (npr = 100) 

x-values for SPL plot (in cm) 
real xlspl , x2spl , x3spl 

parameter (xlspl * 8.5,x2spl - 17, x3spl = 52.5) 


bias velocity 
real vbias 
parameter (vbias = 


= -. 01 ) 


Liner Impedance Parameters 

chc = honeycomb x resistive term 

al = linear resistance term, lower porous sheet 

bl = nonlinear resistance term, lower porous sheet 

c?a = sicrma for resistance distribution . 

saxl - sigma x lower, lower sheet sigma for reactance distribution 

till I s igma~x upper) upper sheet sigrta for reactance drstrrbutrcn 

au = linear resistance term, upper porous sheet 

bu = nonlinear resistance term, upper porous sheet 

xclsp = lower sheet linearized reactance term 

xcusp = upper sheet linearized reactance term 

real chc 

parameter (chc=0 . 8) 
real al,bl 

parameter (al = 0 . 0283 , bl = 21 . 8 ) 
real sg,sgxl,sgxu „ ^ 

parameter (sg = 2 . 0 *dely, sgxl=2 . *sg, sgxu=4 .4*sg) 
real au,bu 

parameter (au=0 . 0905 , bu=233 . 0) 
real xclsp 

parameter (xclsp=- 1 . 14 ) 
real xcusp 

parameter (xcusp=-l . 1338) 



r variable rs is .true. 
T 


if this is a 


restart and .false, if not 


■ variable nmx is the time step number to which to march 
6000 

. variable ck is number of time steps between calls to chkval 0 
100 

. mean ic: .0 top-hat inflow (Not used in LINER) 

2 

* variables dx (= dely = dy) ,dt 

0.1250000 2 . OOOOOOOE-02 

. level of sixth-order artificial dissipation vk 
2 . 0000000E-02 

. constant in reflective outflow condition of Poinsot and Lele (subsonic only 
0 . OOOOOOOE+OO 

, nondime nsional reference pressure for SPL calculation 
1 .4100000E-10 

' * harmonic u disturbance amplitude and frequency 
1 2 . 0000001E-03 0.23099946 


) 

* 


* 

* 

% 

% 

* 

% 

% 

* 

4 

4 

4 

4 

4 

4 

4 


program sp lx . f 

generate SPL for gnuplot macro 'goplot' 

include ' input . inc ' 

real xx(-l :nx+2) , yy ( - 1 :ny+l) 


real spl ( 0 : nx+1 , 0 : ny+1 ) , Y 
integer i, j , il / i 2 , i3 , nl , n2 


open (unit-35 , file- ’ sound' , status-' unknown’ ) 
opSS (unit-351 , file- ' splplot ' , status-' unknown ) 


read (3 5,*) nl,n2 
write (6,*) nl,n2 


do 100 j = 0 , n2-l 

read (35 , * ) (spl (i, j ) / 1=0 

continue 


11 = xlspl/dely + .0000001 

12 = x2spl/dely + .0000001 

13 = x3spl/dely + .0000001 


write (351, 98) xlspl , x2spl , x3spl , 

format ('# xl =',f8.5, x2 - ,f8.5, 


x3 = ' , f 


do 200 j = 0 , n2 -1 
y = j*dely 

write (351, 99) y,spl(il,D 
format (4el3 . 5) 
continue 


) , spl (i2 , j ) , spl (i3 , j ) 


close (35) 
close (351) 





TIME-DOMAIN SIMULATION OF 
ACOUSTIC PROPAGATION IN A LINED DUCT 

R.S. Reichert and S. Biringen Department of Aerospace Engineering Sciences 
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Abstract 

An inviscid, spatial time-domain numerical simulation is employed to compute acoustic wave 
propagation in a duct treated with an acoustic liner. The motivation is to assess the effects on 
sound attenuation of bias flow passed through the liner for application to noise suppression in jet 
engine nacelles. Physically, the liner is composed of porous sheets with backing air cavities. The 
mathematical model lumps the sheets’ presence into a continuous empirical source term which 
modifies the right-hand side of the momentum equations. This source term specifies the time- 
domain characteristics of the frequency-domain resistance and reactance of the liner’s component 
sheets. Nonlinear behavior of the liner sheets at high sound pressure levels is included in the 
form of the source term. The source term constants are empirically matched to frequency-domain 
impedance data via a one-dimensional numerical impedance tube simulation. The resulting liner 
model is then incorporated into a two-dimensional Euler solver and used for simulations of a realistic 
duct configuration. Sound pressure levels and axially transmitted power are computed to assess the 
attenuation effects of various magnitudes of bias flow. Simulation results are compared to available 
experimental data from a geometrically similar lined duct. 

Problem Introduction and Description 

Reduction of noise emitted from jet engines continues to be a key component of aircraft design. 
One major facet is inlet noise. Mechanical and hydrodynamic noise from the engine components 
propagates upstream and out of the inlet. Current generation jet engine designs treat nacelle inlets 
with acoustic liners composed of porous sheets with backing air cavities for acoustic attenuation. It 
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is well known that the impedance of porous sheets varies at high sound pressure levels. Observations 
show that incident sound pressure amplitude relates nonlinearly to particle velocity both for a single 
orifice (Ingard & Ising 1 ) and porous sheets (Melling 2 ). One concept for varying the impedance 
properties of liners is to explort this nonlinear behavior by blowing a low level, steady bias flow 
through them. In this concept, the bias flow, which might be generated naturally by redirecting 
the oncoming air stream through the nacelle skin, can be adjusted to tune the liner impedance for 
optimal attenuation over a range of off-design flight conditions. The simulations presented in this 
work consider high amplitude acoustic propagation within a simple duct geometry in an effort to 
model the physics and eventually validate this idea. 

We implement the liner model into the governing (Euler) equations through a source term in 
the interior of the computational domain. This approach enables the representation of complex 
composite liner structures in terms of one set of empirical, adjustable constants. An alternative 
way of representing the liner is advanced by Tam & Auriault, 3 where finite impedance boundary 
conditions axe imposed at the liner surface. However, because this model requires composite liner 
impedance inf ormation at the boundary, it can be difficult to model complex liners. In the present 
liner model, we lump the physical and geometric attributes of the liner into a source term in the 
momentum equations. This source term is frequency-independent and contains several empirical 
constants which are matched to impedance data (obtained in the frequency-domain), for a given 
liner component. Consequently, the time-domain model so obtained produces the same impedance 
behavior as the actual material. In this way, complex lining geometries, composed of multiple sheets 
with various impedances, can be built up with ease, in contrast to boundary condition methods 
which must know the liner’s composite impedance presented to the duct at the face sheet. The 
governing equations are integrated forward in time to capture the evolution of the acoustic field in 
the presence of the liner structure. The primary motivation is to assess and optimize the effect of 
natural bias flow through the liner. 

The geometry considered here approximates that of a lined duct test section constructed at 
Rohr, Inc., as detailed by Yu, Kwan, & Stockham. 4 Depicted in Fig. 1, the section is 19.25 cm 
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in height and 61 cm in length. The numerical simulation considers a two-dimensional planar cut 
bounded by solid rigid walls at top and bottom and by open boundaries, indicated by dotted lines 
in the figure, at left and right. The experiment also possessed upstream and downstream hard 
walled duct extensions with reverberation chamber terminations; the present simulations mclude 
only the test section. The top surface of the duct is acoustically treated with a liner composed of a 
porous face sheet (18% open area) and backing septum (6% open area) sandwiching a honeycomb 
core. Plane waves (i.e., lowest duct mode) are forced at the left boundary and allowed to propagate 
along the duct to the right through a quiescent air medium. Although only plane wave cases are 
reported here, the method allows examination of higher duct modes simply by changing the form 
of the time-dependent forcing applied at the left boundary. 

As implemented in the present study, time-domain modeling of acoustic lining materials has 
several advantages over frequency-domain analysis. First, it provides a convenient means of imple- 
menting complex liner structure. Composite liners of any number of sheets and backing cavities 
may be built up. As mentioned above, only the component impedance of the sheets, rather than 
the composite impedance at the face sheet, need be known. It is also simple to construct liners m 
which impedance varies spatially, which is useful since segmented treatment allows attenuation of 
widely disparate frequencies (Motsinger & Kraft 5 ). Another strength of the current model is that 
it accommodates both linear and nonlinear noise amplitudes and incorporates nonlinear impedance 
of the porous sheets. Finally, time-domain analysis may treat multiple frequencies and acoustic 
modes simultaneously. These many desirable qualities make time-domain analysis attractive for 
computational aeroacoustic problems. 

Development of Time-Domain Model Form 

Our goal is to develop governing equations for a continuum which contains porous material. 
While no universal form exists for the equations governing flow through porous media (Nayfeh, 
Kaiser, & Telionis 6 ), Morse & Ingard 7 provide a widely accepted form. The discussion below 
justifies, to some degree, the form of their momentum conservation equations. Mass and energy 
conservation equations could have analogous modifications. However, numerical impedance tube 


3 


experiments, similar to those discussed below, reveal that the appropriate mass conservation modi- 
fication produces only infini tesimal differences in results. Energy conservation plays only a passive 
role in acoustics, so its modifications are neglected here as well. Consequently, we treat only 
momentum conservation. ) 

We consider the equations for conservation of momentum, which apply for a compressible, 


viscous flow: 


d , , ^ d , , drtj 

at ( ' n * i ) + aij (w “- ) = 


dxj 


( 1 ) 


for i = 1, 2, 3. Here, Tjj is the stress tensor which includes stress due to pressure. Suppose that 
porous acoustic lining material is distributed uniformly in some region of the flow field. The 


material’s presence will cause flow resistance which can be modeled as as normal stress term across 
an infinitesimal fluid volume in each coordinate direction. This term augments the pressure gradient 
as shown in Fig. 2. Consequently, denoting acoustic material resistance as N, we modify the stress 

term to read 

Tii = -(p + i\0^. (2) 


Here, we have ass umed that the porous material imparts no net shear to the fluid. Following 
Zorumski & Parrott, 8 - 9 we write the gradient of N as a time-domain damping term R t d multiplied 


by the local velocity it*, so that 


drjj 

dxj 


~^--R td««- 

(s%i 


(3) 


Additionally, the effective fluid density within the volume is increased by a constant time-domain 
factor X td (typically between 1.5 and 5.0 for acoustic materials 7 ) due to the material’s presence: 


XtdP 


*td > 1- 


(4) 


The density factor accounts for an increase in effective mass as the fluid moves through constrictions, 
as suggested by Morse & Ingard. 7 Substituting Eq. (3) into Eq. (1) and making the replacement 
(4) yields a modified differential momentum equation: 7 


d_ 

dt 


(pUi) + 


1 

*td 


( dp_ 

\dxi 


+ Rid^t 


)• 


( 5 ) 
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We must now substitute model relations for f? t d and X tc j which properly represent the time-domain 
behavior of acoustic lining materials. This is accomplished by using forms which reproduce the 
porous mat erial behavior in the proximity of each component porous sheet, as elaborated in the 
following section. ) 

Mathematical Model 

The full governing equations are the two-dimensional Euler equations, which express conserva- 
tion of mass, momentum, and energy for the inviscid motion of compressible fluids. They are here 
written for a cartesian domain: 

dW dF dG _ = (6 ) 

-dt + ’ 

where 

W = {p pu pv E) T , ( 7 ) 



\ 


f \ 


pu 


pv 

F = 

pu 2 

G = 

puv 


pvu 


pv 2 


^[E+p]uy 




The state equation closes the system: 


( 8 ) 


E - —2—r + xP(u 2 + v 2 )- 

7-1 2 

The above set has, as the dependent variables, the conserved quantities of mass p, x-directed 
momentum pu, y-directed momentum pv, and total energy E, each expressed on a per volume basis. 
The primitive variables axe density p, x-directed velocity u, y-directed velocity v, and pressure p. 
Note that S is the right-hand side term of Eq. (5). Specification of the open sides as nonreflective 
and the solid walls as rigid slip boundaries in Fig. 1 completes the definition of the mathematical 
problem. 


a 


a 


We have nondimens ionalized the equation system using the following reference scales: 

l T = 0.01 m — > length scale, 

U r = 340.25 m/s — >• velocity scale, 

Ir/Ur = 2.939 x 10 -5 s — > time scale, 
fir = 1.225 kg/m 3 — > density scale, 

p r U r 2 1-418 x 10 5 Pa — >■ pressure, energy 

scale. 

Note that the density and velocity reference values correspond to standard atmosphere density and 
sound speed. It is also worth noting that the equations will capture both linear and nonlinear 
phenomena, which is important for high amplitude wave dynamics. 

The presence of the discrete porous sheets is felt via the source term S which comes from the 
right hand side of Eq. (5). We form R id and X td to capture the general experimental freqwm y- 
domain behavior observed by Ingard k Ising, 1 Melling, 2 and Zorumski k Parrott. 9 Specifically, 
Rohr, Inc. provided impedance data in the form: 

z = {Ho + SV)+ iX, (10) 


where V is the root-mean-square particle velocity at the sheet surface. The real part is the velocity 
dependent resistance, and the imaginary part is the reactance. Reactance is frequency dependent 
and was provided as 


X = mk , 


where k is the incident wavenumber and m is a given constant for each porous sheet type. The 
term 5 (the right- han d side of Eq. (5)) allows for simple nonlinear behavior in the resistance. For 


iitd. we use 


Rtd = fR(x,y){A + B\ui\), 


where Ui = u, v for the x and y momentum equations, respectively. The function /r{x, y) is a 
smearing function meant to distribute the discrete effect of a porous sheet. It is introduced purely 
for numerical convenience and is detailed below. Similarly, we substitute a model for X td in S: 


X td = l-X(f x (x,y)/[f x (x,y)] max)* 
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Here again, we have introduced a convenient smearing function, this time /x(x,y). Once i? t d and 
X t d are specified, we can detail the entire right-hand side source term S: 


V 

o ^ 


0 ) 

1 

[dp/ dx + ufi ?(x, y){A + £|u|)] 


u C(x, y) 

X{fx{x,y)/[fx(x,y)] max) 

[dp/dy -1- vf R (x, y)(A + £|v|)] 


0 


o J 


0 J 


The velocities u and v axe the same particle velocities present in the terms W , F, and G defined by 
Eqs. (7) and (8) and contain the bias flow velocity. In effect, the empirical parameters A, B, and 
X, respectively, specify the levels of linear resistance, nonlinear resistance, and linear reactance of 
ea c h sheet, but in the time-domain. The distribution C[x,y) is set to a constant value (C = 0.8) 
between the face sheet and parallel backing septum to account for the presence of the honeycomb 
core. Since this additional resistance-like term only appears in the x-momentum equation, it acts to 
suppress u within the core and thus helps the face sheet/septum combination approximate a locally 
reacting finer element. Note that A = B = C = X = 0 gives back the original Euler equations, so 
these correspond to open air. 

As discussed above, the functions /j?(x,y) and /x(x,y) control the spatial distribution of the 
liner’s effect for the resistance and reactance, respectively. They are specified as Gaussian distn- 
butions in a coordinate s normal to the sheet surfaces: 

/(.) = I (12) 

Within these functions, a determines the distribution thickness, so we have gr and ax , respectively, 
for Jr{x, y) and /*(x, y). The constant gr is fixed as 2 Ax simply to distribute the discrete effect of 
sheet resistance for numerical purposes. The constants A, B, X, and ax are matched to complex 
impedance data via a numerical impedance tube simulation, as described below. The authors 
conducting the duct experiment 4 provided the impedance data. The above model has a simple 
nonlinear behavior in analogy to Eq. (10), but there is no inherent limitation to the complexity' 
that could be mimicked with the time-domain model. 
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Nonlinear behavior in reactance is not explicitly modeled in Eq. (11), though it is physically 
present. The model consequently linearizes the reactance about each particular bias flow level. 
Nonlinearity tends to reduce the sheet reactance, but the effect is not indefinite as acoustic am- 
plitude increases. When the root-mean-square normal velocity at the sheet surface attains a high 
level, the reactance assumes a constant value. The velocity at which this saturation occurs depends 
on the particular sheet porosity, but is at most V = 0.0095 for the component sheets considered 
here. A constant reactance must consequently be used for all sheets when |vbiasl > 0.0095. Below 
this level, each sheet is considered separately as to whether its reactance behaves nonlinearly or 
is saturated. This means that X in Eq. (10) is a somewhat complicated function of V. To ap- 
proximate the variable reactance, X is consequently computed to be a fixed number between the 
linear X max and the saturated A min , sized according to the expected normal velocity at each sheet 

surface. 

Numerical Solutio n Method 

We discretize and numerically integrate the partial differential equation system Eq. (6), coupled 
with the spatial boundary conditions, in a rectangular domain, with a uniform mesh in both the 
x- and y-directions. The time-advancement scheme is Gottlieb & Turkel’s 10 explicit (2,4) scheme, 
which is second-order accurate in time and fourth-order in space. The variant of the method used 
here is a MacCormack-like explicit predictor/corrector. For grids that are uniform m x and y, the 
method may be written as: 

- FTj) ~ (P&2J - ftijM 

— A y [7(G” J+1 - G"j) - (G?j+ 2 ~ ^ij+i)] 

+At(D"j) + At(^j), ( 13 ^ 

wr / 1 = *{*«+#« 

-Ax[7 (F’j - FUj) - (F*_ij - F:_ 2 j)} 

-A y [7(G* J - - (<%,-_ i - G-j _ 2 )1 
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Here, 


+At(D^) + Af(^)}. 


(14) 


) 



Xy = 


At 

6Ay 


The spatial derivatives of pressure in S are computed using second-order central finite differences. 
The above steps employ forward differencing in the predictor step and backward differencing in the 
corrector step. This is switched on alternating time steps to a predictor with backward differencing 


and corrector with forward differencing to obtain the full (2,4) accuracy. For stability, the Courant- 
Friedrichs-Lewy (CFL) number is made to satisfy: 


CFLx = 


(u + u)max At 2 

Ax 3’ 


and 


CFLy 


[v 4- a)max At ^ 2 

Ay 3’ 


where a is the local sound speed, a = •/ ypj p- 


The (2,4) stencil has five points in both x- and y-directions, so it extends over the domain 
edges at a boundary and first interior point. The method uses the unmodified (2,4) scheme at the 
first interior point. To account for the point beyond the edge, fluxes are extrapolated (third-order) 
to a ghost point one increment outside the domain. Boundary points, including both rigid top 
and bottom walls and the nonreflective boundaries at the “inflow” and “outflow” , are treated with 
the method of Thompson. 11,12 This boundary treatment uses a characteristic decomposition of the 
Euler equations to give an estimate of the flux derivative normal to the boundary at the boundary. 
The method works well with the plane wave modes considered in the present problem, as elaborated 
in the code validation presented below. 

Artificial viscosity terms are added to the numerical scheme to enhance stability, as indicated 
by the D terms. Sixth-order dissipation is added as a source term to the right-hand side of the 
Euler equations: 


D — t<i 



d*W 
d: r 6 


+ (Ay) 6 



(15) 
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The magnitudes of the artificial viscosity terms are 0[(Ax) 6 ] and 0[(Ay) 6 ], so the fourth-order 
spatial accuracy of the method is un a ff ected. The parameter is empirically defined and did not 
exceed = 0.04. The central difference coefficients used for the sixth-derivatives are computed 
using Lagrange polyno m i a l s , found some explicit dissipation necessary to suppress spurious 
oscillations developing around the porous sheets. 

The mean bias flow is specified as a “frozen” field upon which the acoustic perturbations are 
allowed to evolve. In this way, the mean field need not be computed; nor does it evolve. The 


addition of a forcing term e to Eq. (6) accomplishes this: 


dW dF dG * 

~dt + 9* + fy~ S + £ = 0 ' 


(15) 


where 


dF 8G -s 


(17) 


The overbar here denotes the base field, which consists of the bias flow and grazing flow (if present) 
but not the acoustic perturbations. A uniform vertical velocity Ubias is used throughout the domain 


for the base field. Thus, the bias is directed normal to the liner, and there is no grazing flow. For 
the relatively low bias flow levels of the present study, there is physically no signi fi c an t coupling 
between acoustic and mean flow fields, so this frozen flow technique is expected to be valid. It 
should also be noted that a bulk grazing flow in the duct could be implemented using such a 
frozen flow formulation. The technique has the limitation that the mean flow field does not feel the 
presence of the liner structure, but it is a quick way to obtain a converged bias flow. 


Numerical Impedance Tube Simulation 


As mentioned above, the liner model’s parameters must be matched to the a priori known 
impedance data for the component porous sheets under consideration. To t his end, a n um erical 
model for an impedance tube is developed. In this model, one-dimensional Euler equations are 
integrated on a domain with a rigid termination at the right end, and the left boundary is forced 
with acoustic waves. A numerical porous sheet “sample” is placed some distance to the left of 
the rigid termination by centering the Gaussian distributions Eq. (12) at that point. The effective 


impedance of the sample is deduced from the standing wave pattern set up by the incident ;n- J 
reflected waves. 

The impedance tube computations employ parameters similar to those used in the full two- 
dimensional simulations. For instance, Ax = 0.125 and At = 0.025. The acoustic forcing at the 
left boundary of the form: 


p(t) = P cos(ut + n/2) u = p 5 = 0 p = p, 


(18) 


gives purely plane waves with nondimensional wavelength 25, which is also the temporal period. 
This corresponds to 1361 Hz, a value within the frequency range of interest for the two-dimensional 
duct. The domain has 601 points so that the tube length is 75. The simulations are run 5500At, 
or to nondimensional time 137.5. This allows waves propagating at unity sound speed to set up 
several wavelengths of a standing wave without the reflected wave impacting the left boundary. 

We follow the method of Kinsler, Frey, Coppens, & Sanders 13 to extract the impedance from 
the standing wave pattern. Accordingly, the impedance z is computed as 

1 + Ke ie 


z = 


(19) 


1-Ke ie ' 

The constant K is given in terms of the standing wave ratio SWR, the ratio of ma ximum to 
minimum pressure amplitude: 

euro _ 1 D 

( 20 ) 


K . where SWR = Pmax 


SWR+1’ P min 

The phase of the standing waves relative to the sample surface determines 9: 


6 — 2fc(x sam pie ^node) 


( 21 ) 


where k is the waven umb er of the incident waves. (Refer to Fig. 3 for an example of these quantities 
as measured from a standing wave pattern.) The impedance value thus measured equals that of 
the sample plus the impedance of the closed tube behind the sample. Kinsler et al. 13 show that 
this backin g tube impedance is — cot {kd), where d is the depth of the cavity from the sample to 
the rigid ter min ation. Consequently, we fix the depth at one-quarter incident wavelength such that 
— cot(fcd) = 0, and the measured impedance is simply that of the sample. 




A 

B 



crx 

18% 


21.8 


-0.570 

4.0 Ax 

6% 


233 


-0.669 

8.8 Ax 


Table 1: Time-Domain Model Parameters for Liner’s Component Porous Sheets 

* 

Equation (10) represents the frequency-domain impedance, and components 7£o, and X sure 
experimentally determined for a given liner material such that they are known input for matching 
to our time-domain model. The time-domain constants A, B , X , and ax axe adjusted until the 
computed impedance matches that of the porous sheet, given in the form of Eq. (10). In matching, 
there exists a one-to-one correspondence between A and 7£o, B and 5, and X and m = X/k. 
The parameter ax may also be adjusted to effect gross changes in m = X/k. It should be noted 
that X and ax, since they match to m for any wavenumber fc, axe usable time-domain constants 
regardless of the incident frequency. When computing total (linear plus nonlinear) resistance, V is 
computed in the simulation as the root-mean-square velocity over one wave period at the center 
of the sheet distribution. Figure 3 shows a representative standing wave solution for the 6% open 
sheet, with the center of the sheet indicated by the vertical line. Table 1 displays the constants 
obtained through application of this method to the three sheet types. This method is general in 
the sense that complex impedance data for any sample, porous sheet or other, can be matched to 
the empirical constants of the model. 

By matching the frequency-domain data to time-domain parameters, we have ensured that the 
numerical “sample” possesses the correct time-domain behavior. Figure 4 provides an additional 
qualitative validation. The plot shows the Fourier transform of a pressure time series recorded at a 
point midway between the sample and the rigid termination in the impedance tube. This particular 
case employs parameters for a 6% porous sample. It is apparent that the sound pressure transmitted 
through the sample contains overtones of the fundamental almost entirely of odd order. This effect 
is characteristic of material with nonlinear impedance and has been observed experimentally by 
Ingard fe Ising 1 for nonlinear transmission through a single orifice. The numerical liner terms are 
scattering acoustic energy to odd harmonics and sure thus mimicking the correct physical behavior. 
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Duct Code Validation 


The results of two trial simulations provide important checks of the duct code solution method 
and boundary conditions. The solution code has also been extensively tested for Kelvin-Helmholtz 
linpax instabilities down to marginal resolutions, as discussed by Reichert 3 . Those simulations 
showed only slight phase and amplitude errors at resolutions of 12 to 15 points per wavelength. 
The present cases extend the validation to the case of well-resolved (over one hundred points per 
wavelength) acoustic waves. In the first case, a hard wall boundary is placed at the position of the 
face sheet, and the lowest duct mode (i.e., plane wave) is forced with nondimensional amplitude 
0.002 (140 dB ) and frequency 27r/17 (2000 Hz) at the left boundary. The numerical parameters are 
the same as those used for the production duct simulations and are given in the following section. 
In this inviscid s imula tion, the waves should propagate along the duct with no attenuation, so 
sound pressure level (SPL), computed from the simulation data as 


SPL = 10 log 



( 22 ) 


should r emain a constant 140 dB. Here, T is a period of the wave and p Te f is the nondimensional 
equivalent of 2 x 10 -5 Pa. Further, the i-directed intensity I x , defined by 


h = fj 0 “P dt ’ ( 2 3 ) 

should be a nondimensional constant 2 x 10~ 6 at all points within the domain. Finally, the power 
transmitted along the duct per span, computed as 


Pj span = J I x dy , (24) 

should remain constant at 3.05 x 10 -5 nondimensionally, since y max = 15.25. It should be noted 
that, for the calculations of the next section, a dB drop across the duct length is also computed. 
For plane waves (which all later results nearly approximate), it can be shown that the dB drop 
between two x-stations is: 

SPL 2 -SPL! = 10 log/p, (25) 


13 



where f P is the fraction of station 1 power retained at station 2. Plots of the above quantities 
computed from the simulation are not displayed since their profiles are flat. Examination of the 
numerical results, however, shows that the time-domain simulation yields the above quantities to 
at least four significant figures' after reaching a harmonic state. This suggests that the solution 
method and boundary conditions are working properly for this simplified problem. 

The second trial case tests the method’s ability to allow waves to exit the open right boundary 
with minimal reflection. In this case, all significant features of the production simulations are 
considered, including the full duct with a finer composed of a face sheet, septum, and honeycomb. 
Also, a bias flow u bia3 = -0.01 is applied (negative implies blowing out of the liner into the mam 
duct). The duct is forced at the left boundary with 1250 Hz waves of 140 dB. In the first run, the 
domain length is 61, and in the second run, the domain length is doubled. The simulation time 
allows the waves to establish a harmonic state for x < 61 but does not permit them to reach the 
right boundary of the longer domain. This time is also long enough for reflections from the right 
boundary to progress back into the shorter domain. Figure 5 plots transmitted power drop per 
span (from the left side of the domain to the particular x-station of interest) and demonstrates 
that the reflections are negligible and do not degrade the solution in the interior. It is apparent 
that some difference between the curves exists near the outflow, though this difference is less than 
3%. In this study, we use the slope of these curves in their linear region (measured arbitrarily as 
the slope drawn between points at x = 18 and x = 36) as one measure of duct attenuation. The 
two slopes are -2.295 x 10 -7 and -2.297 x 1CT 7 , which differ by much less than 1%. It can be 
concluded that, for the purposes of these simulations, the artificial open boundary conditions are 

performing adequately. 

Duct. Simulation Results 

The effect of varying levels of bias flow is examined for the propagation of the lowest duct 
mode, Eq. (18), within the fined duct of Fig. 1. The domain is discretized using a 489 x 155 
uniform mesh (Ax = Ay = 0.125), and CFL numbers are approximately 0.16 (At = 0.02), which 
is about one-fourth of the numerical method’s linear stability limit. The time-domain solution 
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fields axe marched 6000 At (about two duct length acoustic propagation times), at which time SPL 
anH J x are computed over one harmonic period. The nondimensional amplitude of the waves is 
set to 0.002, which corresponds to 140 dB. The frequency of the waves is u/ = 27r/27.2, which 
dimensionally is 1250 Hz, while bias flow is varied. Nondimensional bias flow velocities are in the 
range |u b iasl < 0.03, which dimensionally gives |vbias| < 10 m/s. Note that this velocity is blown 
out of the liner into the main duct, so it is actually negative in value in the simulation. Other cases 
examine attenuation as acoustic frequency varies over a range from 630 Hz to 4000 Hz at a single 
bias flow level (l^biasl = 0.005). 

Figure 6 is a comparison of sound pressure levels for varying magnitudes of bias flow. The 
plots show SPL versus y at two x-stations. The dotted vertical lines indicate the y-locations of 
the two horizontal porous sheets. As expected, the SPL drops dramatically within the liner as 
l^biasl increases, especially behind the 6% septum. It is apparent that as the bias flow increases, the 
SPL drops most dramatically at the surface of the 6% septum sheet. This is not surprising since 
increasing the bias flow up to |ubiasl = 0.02 changes the resistance of this sheet from its linear value 
of about O.lpc to over 5 pc, which nearly closes off the backing cavity. Even the face sheet attains 
fairly hi gh resistance (about 0.5pc) at the highest bias flow levels depicted in Fig. 6, and the SPL 
shows a (somewhat smaller) drop at its location also. 

Another interesting feature of the SPL plots is that they flatten with increasing bias flow. With 
no bias flow, the SPL tends to show a broad peak near the lower wall (y = 0) and generally 
decreases with increasing y so that it is fairly low near the face sheet (y = 15.25). The application 
of bias flow generally pulls the maximum down at y = 0 and lifts SPL near the face sheet. The 
net effect is a smaller difference between maximum and mi nim u m SPL, or a general flattening of 
the profile. Again, this is consistent with increasing resistance at the face sheet. As discussed in 
the validation section above, a hard wall placed at the face sheet, which has infinite resistance, 
produces completely flat profiles. Thus, it is no surprise to see flatter profiles when Vbias gives the 
face sheet a large resistance. 

Recalling that the hard wall case provides no attenuation of transmitted power suggests that 
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Vbias 

Experimental 
Insertion Loss 
(dB) 

Simulation 

Attenuation 

(dB) 

Experimental 

fp 

% 

Simulation 

fp 

% 

0.0000 

6.125 

5.9 

24.4 

25.5 1 

0.0025 

8.375 

9.1 

14.5 

12.4 

0.0050 

9175 

12.8 

10.6 

5.2 

0.0075 

10.50 

18.1 

8.9 

1.5 

0.0100 

- 11.125 

- 

7.7 

0.0 


Table 2: Comparison of Insertion Loss and Corresponding Fractional Power at Outflow between 
the Present Simulation and the Experiment. 

very high bias flow levels, which increase resistance and make the face sheet act more and more 
like a hard wall, will decrease the liner’s effectiveness. Such a result is seen in Fig. 7. Plotted 
are transmitted power drops for four cases with increasing bias flow magnitudes. It is seen that 
increasing v bias magnitude from 0.0 to -0.01 provides a fairly large decrease in the slope. Over the 
length of this duct, this provides about 40% greater drop in transmitted power. However, increasing 
magnitude from -0.01 to -0.03 provides no additional attenuation. The figure suggests that an 
optimum bias flow level exists for acoustic power attenuation in this lined duct configuration. 

Figure 8 quantifies the optimum bias flow level more clearly. Displayed are the slopes of the 
transmitted power curves, in their linear region, versus bias flow magnitude. Lower (more negative) 
slopes indicate greater power attenuation. It appears that a bias flow tw = “0-01 (3.4 m/s ) pro- 
duces nearly optimum attenuation of transmitted acoustic power down the duct at this frequency. 

Larger bias flows result in a slow decrease in attenuation. 

As mentioned above, the geometry of this study matches the experiment of Yu et al., 4 so that 
the numerical and experimental results may be compared. The experiment measured dB insertion 
loss of a broadband noise source across the duct length at several single frequencies and bias flow 
levels. The maximum attainable bias flow magnitude in the experiment was Kiasl = 0.01. Table 2 
compares attenuations measured at 1250 Hz in the experiment, and their corresponding fraction 
f P of transmitted power remaining by outflow, with the same quantities computed in the present 
study. Both the experiment and present simulation exhibit the same trend of dramatic attenuation, 
followed by a decreasing return, with increasing ]ubiasl- 
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Another comparison with the experiment is depicted in Fig. 9. Values of the simulation’s 
attenuation over the duct length (computed according to Eq. (25)), for Ubias = -0-005, are shown 
at several plane wave frequencies (600 Hz < f < 4000 Hz) in the top plot, while the bottom plot 
is insertion loss as measured in \he experiment. The same frequency response trends are evident in 
both plots: peak attenuation is present at about 1250 Hz with fairly dramatic fall off to either side. 
Note that, on a linear scale, differences when attenuation is 10 dB or more are quite small. The 
approximately 3 dB difference in peak attenuation (at 1250 Hz) corresponds to an experimentally 
observed 95% and a numerically calculated 90% power attenuation. Again, the simulation agrees 
favorably with experiment. 

Concluding Comments 

This study has developed a promising new method for predicting noise attenuation in acousti- 
cally lined ducts. The method is founded upon time-domain governing equations which are modified 
to account for the presence of acoustic lining materials interior to the domain. The modification is 
simply a source term in the momentum equations which provides for acoustic reactance and both 
linear and nonlinear acoustic resistance, but in the time-domain. The modified equations may be 
readily time integrated using numerical methods. The main benefit of analyzing duct problems in 
this way is that complicated liner behavior and structure are easily implemented. Also, multiple 
frequency and multiple mode noise environments may be examined. In sum, the method developed 
here provides promise as an analysis technique for actual liner structures within realistic noise 

environments. 

The time-domain analysis is empirical in the sense that constants in the modified governing 
equations must be matched to complex impedance data for the component sheets of the liner 
using a numerical impedance tube. The impedance tube simulations employ the same equation 
modifications, and the model constants are adjusted until the standing wave pattern in the tube 
matches that which would be produced by a sample with the desired complex impedance. A 
one-to-one correspondence exists between the frequency-domain impedance values and the time- 
domain model constants, which aids in the matching process. Spectral analysis of the pressure 
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signal transmitted through the impedance tube sample gives the same scattering to odd harmonics 
seen in experiment. The matched constants are then used in the full duct simulations to mimic the 

behavior of that material sample. 

Two-dimensional numerical'simulations involving a lined duct provide additional evidence that 
the liner model is a viable-design tool. The computational geometry matches that of an experiment 
conducted at Rohr, Inc. The trends observed in the simulations compare well with those of the 
experiment. Specifically, attenuation as a function of bias magnitude at a single frequency, and 
attenuation as a function of frequency at a single bias magnitude, behave the same in the experiment 
and simulation. The quantitative agreement is acceptable in light of experimental uncertainty, the 
numerical model’s empiricism, and the numerical restriction to two dimensions. It is found that, for 
this liner and 140 dB waves of 1250 Hz , an optimal bias flow magnitude is about 0.01, or 3.4 m/s. 
With the optimum bias velocity, the finer efficiency is significantly improved over the no bias finer, 
lending support to the idea that bias flow may be used to tune and improve the performance of a 
liner. Further, this application demonstrates that the new model may be used as a design tool for 
determining optimal finer configurations and operating conditions. 
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Figure 3: Numerical Impedance Tube Solution 


Pressure Amplitude 



Figure 4: Spectrum of Pressure Time Series between Sample and Rigid Termination in Impedance Tube 
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Figure 9: Attenuation Frequency Response Comparison: Simulation and Experiment 
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Abstract 

An invisdd, spatial time-domain numerical simu- 
lation is employed to compute acoustic wave prop- 
agation in a duct treated with an acoustic liner. 
The motivation is to assess the effects on sound at- 
tenuation of bias flow passed through the liner for 
application to noise suppression in jet engine na- 
celles. Physically, the liner is composed of porous 
sheets with backing air cavities. The mathemat- 
ical model lumps the sheets’ presence into a con- 
tinuous empirical source term which modifies the 
right-hand side of the momentum equations. This 
source term specifies the time-domain behavior of 
the frequency-domain resistance and reactance of 
the liner’s component sheets. The source term con- 
stants are matched to frequency-domain unpedance 
data via a one-dimensional numerical impedance 
tube simulat ion. Nonlinear behavior of the liner at 
high sound pressure levels is included in the form 
of the source term. Sound pressure levels and axi- 
ally tr ansmi tted power are computed to assess the 
effect of various magnitudes of bias flow on attenu- 
ation. Simulation results are compared to available 
experimental data on a geometrically similar lined 
duct. 

Problem Introduction and Description 

Reduction of noise emitted from jet engines con- 
tinues to be a key component of aircraft design. One 
major facet is inlet noise. Mechanical and hydro- 
dynamic noise from the engine components propa- 
gates upstr eam and out of the inlet. Current gen- 
eration jet engine designs treat nacelle inlets with 
acoustic liners composed of porous sheets with back- 
ing air cavities for acoustic attenuation. It is well 
known that the impedance of porous sheets var.es 
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at high sound pressure levels. Observations show 
that incident sound pressure amplitude relates non- 
lineaxly to particle velocity both for a^ single orifice 
(Ingard 1 ) and porous sheets (Melling 2 ). One con- 
cept for varying the impedance properties of liners is 
to exploit this nonlinear behavior by blowing a low 
level, steady bias flow through them. In this con- 
cept,' the bias flow, which might be generated nat- 
urally by redirecting the oncoming stream through 
the nacelle skin, could be adjusted to tune the liner 
impedance for optimum attenuation over a range of 
flight conditions. The present study considers high 
amplitude acoustic propagation within a simple duct 
geometry in an effort to model the physics and even- 
tually validate this idea. 

The goal of this work is to develop a design tool 
for time-domain numerical simulation of acoustics 
in the presence of sound-absorbing liners with bias 
flow. The physical and geometric attributes of the 
liners are lumped into a model term with empiri- 
cal constants. This model term modifies the gov- 
erning equations such that the liner manifests itself 
interior to the domain rather than through finite 
impedance boundary conditions. The model con- 
stants are matched to complex impedance data for 
liner materials so that the simulation produces the 
same impedance behavior as the actual material. 
The governing equations are integrated forward in 
time to capture the evolution of the acoustic field in 
the presence of the finer structure. The primary mo- 
tivation is to assess and optimize the effect of natural 
bias flow through the finer. 

The geometry considered here approximates that 
of alined duct test section constructed at Rohr, Inc., 
as detailed in Yu, Kwan, k Stockham 3 . Depicted 
in Fig. 1, the section is 19.25 cm in height and 
61 cm in length. The numerical simulation consid- 
ers a two-dimensional planar cut bounded by solid 
rigid walls at top and bottom and by open bound- 
aries, indicated by dotted lines in the figure, at left 
and right. The top surface of the duct is acousti- 
cally treated with a liner composed of a porous face 
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sheet (18% open area) and backing septum (6% open 
area) san dwiching a honeycomb core. Additionally, 
the air cavity behind the septum may contain back- 
ing sheets (8.7% open area). These sheets approx- 
imate the presence of corrugation septa which are 
distributed in the duct’s third dimension in the ac- 
tual experiment; two configurations distributing the 
septa in the simulation’s two dimen s i ons were tested, 
as will be discussed. Plane waves (i.e., lowest duct 
modes) are forced at the left boundary and allowed 
to propagate along the duct to the right through a 
quiescent air medium. 

Development of Time-Domain Model Form 

Consider the integral form of in viscid conservation 
of momentum applied to a fixed control volume, 


[ ?k?LdV = - [{pv-dA)v- f pdA, 

Jv dt J a •>* 


( 1 ) 


For the differential control volume- 
flux term becomes: 

a 


momentum 


which states that the rate of change of momentum 
within the control volume is equal to the net flow of 
momentum across the volume’s surface plus the mo- 
mentum change due to surface pressure forces. Sup- 
pose now that we apply this equation to the Carte- 
sian differential control volume = dV = 

dxdydz) depicted in Fig. 2. Further, suppose that 
some porous acoustic linin g material is distributed 
uniformly throughout the volume. The material’s 
presence will lead to a pressure change across the 
volume in osch coordinate direction, augmenting the 
pressure gradient as shown. The pressure force term 
becomes 

<2) 

where we have employed tensor notation on the right 
side. Following Zorumski Parrott 4, , we write this 
pressure jump term as a time-domain resistive term 
Rt d multiplying the local velocity u*, 


(£).- 




so that 


lim —Ip dA — 
AV-.0 J A 


_|£. dV - fltdUi dV. 

dxi 


lim — f (jn-dA) « = *td-£~ >=* v i) . (6) 

a v — *o J A CrZ- 


hmo Eq. (1) and 
modified differ- 


Substituting Eqs. (4), (5), and (6) 
recovering vector notation yields a 
ential momentum equation: 

^l + V-(ptf5) = ~(Vtr-Jl td tO. (7) 

Ot 

We mus t now substitute model reLsmms for and 
Xtd which properly mimic the time— a i n behavior 
of acoustic lining materials. 

Mathematical Moie 

The gover nin g equations are the Two-dimensional 
Euler equations, which express conservation of mass, 
momentum, and energy for the mviscid motion of 
compressible fluids. They axe here written for a 
Cartesian domain with no body forces present: 

( 8 ) 


(3) 


(4) 
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at + dx dy 
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( 9 ) 


The state equation doses the system: 

P 


E = 


7-1 


4- |p(« 2 + v 2 ). 


( 10 ) 


( 11 ) 


Equation (3) assumes that the pressure jump will 
instantaneously satisfy a steady flow resistance rela- 
tion across the differential volume. Additionally, the 
effective fluid density within the volume is reduced 
by a constant time-domain factor Xui between 0 and 
1 due to the material presence: 

p — ► X^p 0 < X td < 1- (5) 


The above set has, as the dependent variables, the 
conserved quantities of mass p, x-directed momen- 
tum pu, y-directed momentum pv, and total energy 
£, each expressed on a per volume basis. The prim- 
itive variables are density p, x-directed velocity u, 
y-directed velocity u, and pressure p. Specification 
of the open sides as nonreflective and the solid walls 
as rigid slip boundaries completes the definition of 
the mathematical problem. 

We have nondimensionaiized the equation system 
using the following reference scales: 


l r = 0.01 m 
U r - 340.25 m/s 
IrfUr = 2.939 X 10“ 5 s 
p r = 1.225 kg/m? 
p r Ur 2 = 1.418 X 10 5 Pa 


length scale, 
velocity scale, 
time scale, 
density scale, 
pressure, energy 
scale. 
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Note that the density and velocity reference val- 
ues correspond to standard atmosphere density and 
sound speed. It is also worth noting that the equa- 
tions will capture both linear and nonlinear phenom- 
ena. 

The presence of the discrete porous sheets is felt 
via the source term S from the right hand side of 
Eq. (7). We form Rtd and X t a to capture the general 
experimental frequency-domain behavior observed 
by Ingard 1 , Melling 3 , and Zo ruins ki & Parrott 5 . 
Specifically, Rohr, Inc. provided impedance data in 
the form 

z = (TZq -I- 5V) *f t X , (12) 

where V is the root-mean-square velocity at the 
sheet surface. The real part is the velocity depen- 
dent resistance, and the imaginary part is the re- 
actance. The term S allows for simple nonlinear 
behavior in the resistance: 

- 1 

S ~ ~l-X{fx(.x,y)/[fx(x,y)] max) 


0 \ 

[dp/dx + u/r(x , y)(A + S|u|)] 

[dp/dy + vfn(x , y)(A -h B\v\)] 

0 / 


f 0 \ 

t! C{x,y) 

0 

V 0 / 


(13) 


The em pirical parameters A, B, and X respectively 
specify the levels of linear resistance, nonlinear re- 
sistance, and linear reactance of each sheet, but in 
the time-domain. The distribution C(x,y) is set 
to a constan t value (C = 0.8) between the free 
sheet and parallel backing septum to account for 
the presence of the honeycomb core. Since this ad- 
ditional resistance-like term only appears in the x- 
momentum equation, it acts to suppress it within 
the core and thus helps the face sheet/septum com- 
bination appro ximat e a locally reacting liner surface. 
Note that A=B=C=X= 0 gives back the orig- 
inal Euler equations, so these correspond to open 
air. The functions /ji(x,y) and /x(x,y) control the 
spatial distribution of the liner’s effect for the resis- 
tance and reactance, respectively. They are specified 
as Gaussian distributions in a coordinate s normal 
to the sheet surfaces: 




(14) 


The constant <tr is fixed as 2Ax simply to distribute 
the discrete effect of sheet resistance for numeri- 
cal purposes. The constants .4, B , X , and crx axe 
matched to complex impedance data via a numeri- 
cal impedance tube simula tion, as described below. 
The authors conducting the duct experiment 3 pro- 
vided the impedance data. 

Nonlinear behavior in reactance is not explicitly 
modeled in Eq. (13). This nonlinearity tends to re- 
duce the sheet reactance, but the effect is not in- 
definite in the provided impedance data. When the 
root-mean-square normal velocity at the sheet sur- 
face attains a high enough level, the reactance again 
assum es a constant value. The velocity at which this 
saturation occurs depends on the particular sheet 
porosity, but is at most 0.0095 for the component 
sheets considered here. A constant reactance must 
be used for all sheets when |vbias| ^ 0.0095. Below 
this level, each sheet must be considered separately 
as to whether its reactance behaves nonlinearly or is 
saturated. 

It was found in the course of the simulations that 
if an explicit nonlinear reactance term is used in 
Eq. (13), the model produces a smaller reactance re- 
duction thAn it should. This resulted in large jumps 
in the averaged data (i.e., in the slope of the trans- 
mitted power per span curves seen in Fig. 7) when 
the bias flow magnitude dictated that the reactance 
switch from nonlinear to saturated levels. Of course, 
this only affected cases for which |vbias| < 0.0095, 
since all of the sheets perform with saturated re- 
actance for (tw| > 0.0095. To approximate the 
variable reactance, X was computed to be a fixed 
number between the X max and the saturated 

Xaunj sized according to the expected normal veloc- 
ity at each sheet surface. Effectively, this process is a 
linearization of the reactance about each particular 
bias flow level. 

Nirnifiriral Solution Method 

We wish to discretize and integrate the partial 
differential equation system (8), coupled with the 
spatial boundary conditions. The do ma in is rect- 
angular, with uniform mesh in both the x- and y- 
directions. The time- advancement scheme is Got- 
tlieb & Tinkers 6 explicit (2,4) scheme, which is 
second-order accurate in tune and fourth-order in 
space. The variant of the method used here is a 
MacCormack-like explicit predictor/corrector. For 
grids that are uniform in x and y, the method may 
be written as: 


529 


*&' 


-A.[7(j S-w - fly - tfllu - ^w)l 
— A Jf (7(G7 J+1 — G£,-) — (<5^,+2 — ^£i+i)] 
+Af(D)^ + Ai(5)r j , (15) 

-A,[7(FV - i$L W ) - (3-W " ^-2,i)l 
-A y [7(GV - GV_ t ) - (GVi - <?Jj-a)] 

+Af(D)*j 4- Af(S)'j (16) 


A s = A v = ^-. 

The above steps employ forward differencing in the 
predictor step and backward differencing in the cor- 
rector step. This is switched on alternating time 
steps to a predictor with backward differencing and 
corrector with forward differencing in order to obtain 
the full (2,4) accuracy. For stability, the Courant- 
Friedrichs-Lewy (CFL) number should satisfy: 


CFL X = 


CFLy — 


(u + a) 


a)maxA t 2 

Ax 3’ 


(v + a) c 


where a is the local sound speed, a = y/jp/p- 

The (2,4) stencil has five points in both x- and 
y-directions, so it extends over the domain edges at 
a boundary and first interior point. The method 
uses the unm o difi ed (2,4) scheme at the first inte- 
rior point. To account for the point beyond the edge, 
fluxes are extrapolated (third-order) to a ghost point 
one increment outside the domain. Boundary points 
are treated with the method of Thompson 7,8 . This 
boundary treatment uses a characteristic decompo- 
sition of the Euler equations to give an estimate of 
the flux derivative normal to the boundary at the 
boundary. 

Artificial viscosity terms may be added to the 
scheme to enhance stability, as indicated by the D 
terms above. Sixth-order dissipation is added as a 
source term to the right-hand side of the Euler equa- 




The sizes of the artificial viscosity terms are 
0((Ax) 6 ] and 0[(Ay) 6 ], so the fourth-order spatial 


accuracy of the method is unaffected. The parame- 
ter td is user-defined and did not exceed td — 0.04. 
The central diff erence coefficients used for the sixth- 
derivatives are computed using Lagrange polynomi- 
als. We found some explicit dissipation necessary to 
suppress spurious oscillations developing around the 
porous sheets. 

The mean bias flow is sperified as a “frozen” field 
upon which the acoustic perturbations are allowed 
to evolve. In this way, the mean field need not be 
computed; nor does it evolve. The addition of a 
forcing term Z to Eq. (8) accomplishes this: 

dW dF , dG 9 - 

dt dx dy 


where 


dF dG ^ 
€ ~ dx dy 


The overbar here denotes the mean field. A uniform 
vertical velocity is used throughout the domain 
for the mean field. Thus, the bias is directed nor- 
mal to the liner, and there is no grazing flow. This 
technique has the limitation that the mean field does 
not feel the presence of the liner structure, but it is 
a quick way to obtain a converged bias flow. 

N"Tneri<-al Tmpedance Tube S imula tion 

As mentioned above, the liner model’s empirical 
parameters must be matched to the impedance data 
for the porous sheets under consideration. To this 
end, a numerical model of an impedance tube wa<> 
developed. One-dimensional Euler equations are in- 
tegrated on a domain with a rigid te rmi n a tion at the 
right end. The left bo undar y is forced with acoustic 
waves. A numerical porous sheet “sample” is placed 
some distan ce to the left of the rigid termination 
by centering the Gaussian distributions (14) at that 
point. The effective impedance of the sam ple is de- 
duced from the standing wave pattern set up by the 
incident and reflected waves. 

De tails of the impedance computation may be 
found in introductory acoustics texts. Here, we fol- 
low that of Kinsler, Ftey, Coppens, ic Sanders 9 . The 
measured impedance z is computed as 

( 20 ) 

1 - Ke‘ a 

The constant K is given in terms of the standing 
wave ratio SWR, the ratio of maximum to minimum 
pressure amplitude: 


SWR - 1 


where SWR = 



The value for 9 is given by the phase of the standing 
waves relative to the sample surface: 

0 = 3?node) C^) 

where k is the wavenumber of the incident waves. 
The impedance value thus measured equals that of 
the sample plus the impedance of the closed tube 
behind the sample. Kinsler et al. 9 show that this 
backing tube impedance is — cot(id), where d is the 
depth of the cavity from the sample to the rigid ter- 
mination. Consequently, we fix the depth at one- 
quarter incident wavelength such that — cot (fed) = 
0, so that the measured impedance is simply that of 
the sample. 

The tune-domain constants A, B 1 X , and ax 
are adjusted until the computed impedance matches 
that of the porous sheet, which is of the form of 
Eq. (12). The present model has a simple nonlinear 
frequency-domain behavior, but there is no inherent 
limitation to the complexity that could be mimicked 
with a time-domain model. In matching, there ex- 
ists a one-to-one correspondence between A and TiLq, 
B and 5, and X and X. The parameter ax may 
also be adjusted to effect gross chan ges in X. Fig- 
ure 3 shows a representative standing wave solution 
for the 6% open sheet, with the center of the sheet 
indicated by the vertical line. The following table 
shows the constants obtained through application of 
this method to the three different sheet types. 



This method is general in the sense that complex 
impedance data for any sample, porous sheet or 
other, can be matched to the empirical constants 
of the model. 

Duct Code Validation 

The results of two trial simulations provide impor- 
tant checks of the duct code solution method and 
boundary conditions. In the first case, a hard wall 
boundary condition is placed at the position of the 
face sheet, and the lowest duct mode is forced with 
nnn dimens ional amplitude 0.002 (140 dB) and fre- 
quency 27r/17 (2000 Hz) at the left boundary. In 
this inviscid simulation, the waves should propagate 
along the duct with no attenuation, so sound pres- 
sure level (SPL), computed from the simulation data 
as 

SPL = 10 l °9\(fJ o P 2 A fa, - 


should remain a constant 140 dB. Here, T is a period 
of the wave and p rc / is the nondimen sional equiva- 
lent of 2 x 10~ 5 Pa. Further, the x-directed intensity 
/*, defined by 

ix= fJ Q dt ' ( 24 ) 


should be a nondimensional constant 2 x 10“ 6 at 
all points within the domain. Finally, the power 
transmitted along the duct per span, computed as 



should remain constant at 3.05 x 10“ 5 nondimcb 
sionally, since ymax — 15.25. Plots of the above 
quantities as computed from the simulation are not 
displayed since all profiles are flat. Examination of 
the numerical results, however, show that the time- 
domain simulation yields the above quantities to at 
least four significant figures after reaching a har- 
monic state. 

The second trial case tests the method’s ability to 
allow waves to exit the open right boundary with 
minimal reflection. In this case, all si gnifi c a n t fea- 
tures of the simulations are included. That is, the 
fall domain is considered, including the presence of a 
liner with face sheet, septum, and honeycomb. Also, 
a bias flow = —0.01 is applied. One run em- 
ploys the domain length of 61 to be used for the 
later simulations; the domain length is doubled for 
the second run. The time of simulation is designed so 
that the waves have a chance to establish a harmonic 
state for x < 61 but do not have a chance to reach 
the right boundary of the long domain. The time 
is also long enough for reflections to appear in the 
shorter domain. Figure 4 demonstrates that degra- 
dation from reflections is slight. Plotted are drops 
in transmitted power per span from the left side of 
the domain to the particular x-station of interest. It 
is apparent that some difference between the curves 
exists near the outflow, though this difference is less 
than 3%. In this study, we use the slope of these 
curves in their linear region (measured arbitrarily 
as the slope drawn between points at x = 22 and 
x = 35) to quantify the duct attenuation. The two 
slopes axe -2.035 x 10“ 7 and 2.034 x 10“ 7 , which 
differ by less than 1%. It can be concluded that for 
the purposes of these simulations, the artificial open 
boundary conditions are performing adequately. 


(23) 


Duct Simulation Results 

The effect of varying levels of bias flow is exam- 
ined for the propagation of the lowest duct mode 









when the liner consists of the face sheet and parallel 
septum sandwi chin g a honeycomb core. The domain 
is discretized using a 489 x 155 uniform mesh (Ax = 
Ay = 0.125), and CFL numbers are approximately 
0.16, which is about one-fourth of the linear stability 
limit. The nondimensional amplitude of the waves 
is set to 0.002, which corresponds to 140 dB .. The 
frequency of the waves is cj = 2*r/17, which dimen- 
sionally is 2000 Hz. Nondimensional bias flow ve- 
locities are in the range 0.001 < |vbiasl - 0.03, which 
dimensionally gives 0.34 mj s < [ubiasl - ^ m/s. 

Figure 5 is a comparison of sound pressure levels 
for varying magnitudes of bias flow. The plots show 
SPL cuts in y at two x-stations. The y-location of 
the two horizontal porous sheets is indicated by the 
dotted vertical lines. As expected, SPL drops dra- 
matically within the liner as |vbi*s| increases. It is 
apparent that as the bias flow increases, SPL drops 
more dr amat ically at the surface of the 6% septum 
sheet. This is sensible because increasing bias flow 
up to [vbiasl = 0.02 changes the resistance of this 
sheet from its linear value of about O.lpc to over 5 pc, 
which essentially closes off the ba cki ng cavity. Even 
the face sheet at tains fairly high resistance (about 
0.5pc) at the highest bias flow levels depicted in Fig. 

5, and SPL shows a (somewhat smaller) drop at its 
location also. 

Another interesting feature of the SPL plots is 
that they flatten with increasing bias flow. With 
no bias flow, SPL tends to show a broad peak near 
the lower wall (y = 0) and generally decreases with 
increasing y so that it is fairly, low at the face sheet 
(y = 15.25). The application of bias flow generally 
pulls the down at y = 0 and lifts SPL 

near the face sheet. Again, this is consistent with 
increasing resistance at the face sheet. As discussed 
in the validation section above, a hard wall placed at 
the face sheet, which is infinite resistance, produces 
completely flat profiles. Thus, it is no surprise to see 
fairly fiat profiles when Vbias gives the face sheet a 
large resistance. 

Recalling that the hard wall case provided no at- 
tenuation of tr ansmi tted power suggests that very 
high bias flow levels, which make the face sheet act 
more and more like a hard wall, will decrease the 
liner's effectiveness. Such a result is seen in Fig. 6. 
Plotted are transmitted power drops for four cases 
with increasing bias flow magnitudes. It is seen that 
increasing Vbu* from 0.0 to -0.01 provides a fairly 
large decrease in the slope. Over the length of this 
duct, this provides about 30% greater drop in trans- 
mitted power. However, increasing tfbias from -0.01 
to -0.02 provides no additional attenuation. The fig- 
ure suggests an optimum bias flow level for acoustic 


power attenuation in this duct /liner configuration. 

Figure 7 q uantifi es the- optimum bias flow level 
more clearly. Displayed are the slopes of the trans- 
mitted power curves in their linear region versus bias 
flow magnitude. Lower slopes indicate greater power 
attenuation. The first plot exhibits slopes for the 
basic configuration of an 18% face sheet plus 6% 
parallel b ackin g septum with an empty backing cav- 
ity. For the second plot, cases were run with 8.^ %, 
±45° corrugation sheets present in the ba ckin g cav- 
ity, but running in the x-direction. Finally, cases for 
the third plot were run with an 8.7% parallel sep- 
tum present in the backing cavity midway between 
the 6% septum and the rigid wall. The latter two 
configurations are attempts to, in some measure, ac- 
count for the presence of 8.7% sheet corrugations 
distributed in the third dimension of the actual ex- 
perimental duct lining 3 . All cases include the hon- 
eycomb between the face sheet and septum. It is 
seen in Fig. 7 that the presence of 8.7% sheets in 
the backing cavity only marginally affects the mini- 
mum slope. For all three configurations, it appears 
that a bias flow t* ias = -0.01 (3.4 m/s) produces 
optimum attenuation of tr ansm itted acoustic power 
down the duct at this frequency. 

As mentioned above, the geometry of this study 
matches the experiment of Yu et al. 3 The exper- 
iment measured dB insertion loss of a broadband 
noise source across the duct length at several differ- 
ent frequencies and bias flow levels. The following 
table compares attenuations measured at 2(300 Hz 
in the experiment to the power slopes computed in 
the present study; 



Experiment 

Simulation 


Insertion Loss 

Power Slope 


(dB) 

(xlO- 7 ) 

0.0000 

3.50 

-1.35 

0.0025 

4.25 

-1.60 

0.0050 

5.25 

-1.70 

0.0075 

5.25 

-1.S0 

0.0100 

4.75 

-1.85 

0.0125 

— 

-1.85 

0.0250 

- 

-1.75 


Both the experiment and present simulations exhibit 
the same trend of rapidly increasing attenuation, 
followed by slowly decreasing attenuation, with in- 
creasing jvbiasl* The opt im u m ut>i*s hi the experi- 
ment has a somewhat lower magnitude, but overall , 
the agreement is favorable. 

Another comparison is depicted in Fig. 8. Values 
of the simulation's power slope for Vb ia5 = —0.005 
are shown at several different plane wave frequen- 
cies in the upper plot, while the lower plot shows 
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insertion loss as measured in the expe rimen t. The 
same frequency response trends are evident in both 
plots: peak attenuation is present at about 1250 Hz 
with fairly dramat ic fall off to either side. Again, 
the simu lation agrees favorably with experiment. 

Con cluding Comments 

This study has developed a new method for pre- 
dicting noise attenuation in acoustically lined ducts. 
Time-domain simulation with governing equations 
modified to include the time-domain behavior of 
liner resistance and reactance was applied to a duct 
test section geometry. A liner configuration com- 
posed of a face sheet and parallel septum surround- 
ing a honeycomb core of depth 1-25 cm was set on 
one side of the 15.25 cm duct. A cavity of depth 
2.75 cm backed the septum; distribution of addi- 
tional septa within this backing cavity produced 
minim al change in the results. The time-domain 
liner model included nonlinear behavior in the re- 
sistance of the liner’s component sheets. The sim- 
ulations have demonstrated the presence of an op- 
timum bias flow level. For waves of 140 dB and 
frequency 2000 Hz propagating axially along the 
duct, |vbias| = 3.4 m/s produced the greatest at- 
tenuation in tr ans mitted acoustic power along the 
duct. Trends with varying bias flow and varying 
frequency agreed favorably with experimental mea- 
surements performed at Rohr, Inc. 
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Figure 3: Numerical Impedance Tube Solution 
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Figure 4: Axial Power per Span Drop from Left Boundary: Short vs. Long Domain. 
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Figure 5: SPL Levels for Increas i ng Bias Flow Magnitudes 
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Figure 6: Axial Power per Span Drop for Increasing Bias Flow Magnitudes 





Slope 





Figure 7: Slope of Axial Power per Span in Linear Region for Three Liner Configurations 







